Graphical display of raster data

Advanced Geospatial Data Processing for Social Scientists

Dennis Abel & Stefan Jünger

2025-04-28

Now

Day Time Title
April 28 10:00-11:15 Introduction
April 28 11:15-11:30 Coffee Break
April 28 11:30-13:00 Raster data in R
April 28 13:00-14:00 Lunch Break
April 28 14:00-15:15 Raster data processing
April 28 15:15-15:30 Coffee Break
April 28 15:30-17:00 Graphical display of raster data in maps
April 29 10:00-11:15 Datacube processing I
April 29 11:15-11:30 Coffee Break
April 29 11:30-13:00 Datacube processing II & API access
April 29 13:00-14:00 Lunch Break
April 29 14:00-15:15 Data integration and linking (with survey data)
April 29 15:15-15:30 Coffee Break
April 29 15:30-17:00 Outlook and open session with own application

What is ggplot2?

In general, ggplot2 is well-suited for multi-dimensional data like raster layers and stacks.

Components of the plot are added as layers.

plot_call +
  layer_1 +
  layer_2 +
  ... +
  layer_n

If you are new to ggplot2, you might want to check out:

Components of a Plot

According to Wickham (2010, p. 81), a layered plot consists of the following components:

  • Data and aesthetic mappings,
  • Geometric objects,
  • Scales,
  • (and facet specification)
plot_call +
  data +
  aesthetics +
  geometries +
  scales +
  facets

Recap: Plotting vector data

Load the data

# load district shapefile
german_districts <- sf::read_sf("./data/VG250_KRS.shp")

# load district attributes
attributes_districts <- 
  readr::read_csv2("./data/attributes_districts.csv") |> 
  dplyr::mutate(ecar_share = as.numeric(ecar_share))

# join data
german_districts_enhanced <- 
  german_districts |>  
  dplyr::left_join(attributes_districts, by = "AGS")

# load states shapefile
german_states <- sf::read_sf("./data/VG250_LAN.shp")

Adding geoms to a blank canvas

# a simple first map 
ggplot(data = german_districts_enhanced)

Adding geoms to a blank canvas

# a simple first map 
ggplot() +
  geom_sf(data = german_districts_enhanced)

Add the aesthetics

Are you having trouble choosing the right color? Some excellent tutorials exist, f.e, by Michael Toth.

# change color palette
ggplot() +
  geom_sf(
    data = german_districts_enhanced, 
    aes(fill = ecar_share)
  ) + 
  # readable with color vision deficiencies
  scale_fill_viridis_c(option = "plasma", direction = -1)

Add another layer

# the shapefile includes polygons of oceans and lakes
# easy fix on the fly when you know your data
german_states <-
  german_states |>  
  dplyr::filter(GF == 4)

# add layer with German states
ggplot() +
  geom_sf(
    data = german_districts_enhanced, 
    aes(fill = ecar_share), 
    color = NA
  ) + 
  scale_fill_viridis_c(
    option = "plasma", 
    direction = -1
  ) +
  # add another layer
  geom_sf(
    data = german_states, 
    # filling transparent
    fill = "transparent",
    # color of borders
    color = "black", 
    # size of borders
    size = 1
  )

Save and reuse

Maps produced with ggplot2 are standard objects like any other object in R (they are lists). We can assign them to reuse, plot later, and add map layers.

Furthermore, you can save them just as any ggplot2 graph. The ggsave() function automatically detects the file format. You can also define the height, width, and dpi, which is particularly useful to produce high-class graphics for publications.

Save and reuse

# assign to object
ecar_map <- 
  ggplot() +
  geom_sf(
    data = german_districts_enhanced, 
    aes(fill = ecar_share), 
    color = NA
  ) + 
  scale_fill_viridis_c(
    option = "plasma",
    direction = -1,
    name = "E-Car Share",
    guide = guide_legend(
      direction= "horizontal",
      label.position = "bottom"
    )
  ) + 
  geom_sf(
    data = german_states, 
    fill = "transparent", 
    color = "black"
  ) 

# save as png-file
# ggsave("ecar_map.png", ecar_map, dpi = 300)

Where ggplot2 cannot help anymore

In some specific circumstances, we might realize that ggplot2 is super powerful but not originally designed to build maps. Typical features of maps are not in the package, like a compass or scale bars.

This is where other packages might need to be installed. The good thing: Elements of the package ggspatial can be included as ggplot2 layer. Check out Github.

The extras

ggspatial allows you to add, f.e. a scale bar and a north arrow.

# add scalebar and north arrow
ecar_map +
  ggspatial::annotation_scale(
    location = "br"
  ) +
  ggspatial::annotation_north_arrow(
    location = "tr", 
    style = ggspatial::north_arrow_minimal()
  )

Making a plan

This map will be our canvas for the ongoing session. There are hundreds of options to change this map. We will cover at least some essential building blocks:

  • THE MAP: adding attributes, choosing from colors/palettes, adding layers
  • THE LEGEND: position, sizes, display
  • THE ENVIRONMENT: choosing from themes and build your own
  • THE META-INFORMATION: titles and sources
  • THE EXTRAS: scales and compass

If you are working on your maps, the ggplot2 cheatsheets will help you with an overview of scales, themes, labels, facets, and more.

ggplot2 and raster data

You can perfectly use ggplot2 to create maps with raster data. There are several ways to do so. The easiest way is using the tidyterra package.

library(tidyterra)

# Create random raster
r <- terra::rast(nrows = 50, 
                 ncols = 50, 
                 xmin = 0, 
                 xmax = 10, 
                 ymin = 0, 
                 ymax = 10
                 )
terra::values(r) <- runif(terra::ncell(r))

# Plot with tidyterra::geom_spatraster
ggplot() +
  geom_spatraster(data = r) +
  scale_fill_viridis_c(na.value = "transparent") +
  theme_minimal()

Case study: Population dynamics

Let’s explore the suitability of ggplot2 in combination with tidyterra with a case study on population dynamics in Germany. We are utilizing the population grids from the WorldPop Open Population Repository (WOPR).

# Example with worldpop data for Germany 2020
ger_pop_2020 <- terra::rast("../../data/deu_ppp_2020_1km_Aggregated.tif")

ger_pop_2020
class       : SpatRaster 
dimensions  : 934, 1100, 1  (nrow, ncol, nlyr)
resolution  : 0.008333333, 0.008333333  (x, y)
extent      : 5.87375, 15.04042, 47.27458, 55.05792  (xmin, xmax, ymin, ymax)
coord. ref. : lon/lat WGS 84 (EPSG:4326) 
source      : deu_ppp_2020_1km_Aggregated.tif 
name        : deu_ppp_2020_1km_Aggregated 
min value   :                       0.000 
max value   :                    4513.521

Simple plot with terra/base R

terra::plot(ger_pop_2020)

Simple plot with terra/base R

Let’s transform it into ETRS89/UTM 32N (EPSG: 25832). You will detect a slight adjustment to the visualization.

ger_pop_2020 <- terra::project(ger_pop_2020, 
                               "EPSG:25832"
                               )

terra::plot(ger_pop_2020)

Let’s turn to ggplot

ggplot()+
  geom_spatraster(data = ger_pop_2020)

Let’s turn to ggplot

Adjust the color scheme

ggplot()+
  geom_spatraster(data = ger_pop_2020)+
  scale_fill_whitebox_c(
    palette = "muted",
    n.breaks = 12,
    guide = guide_legend(reverse = TRUE)
  )

Let’s turn to ggplot

Adjust the color scheme

ggplot()+
  geom_spatraster(data = ger_pop_2020)+
  scale_fill_viridis_c(
    n.breaks = 12,
    guide = guide_legend(reverse = TRUE)
  )

Let’s turn to ggplot

Add labels

ggplot()+
  geom_spatraster(data = ger_pop_2020)+
  scale_fill_viridis_c(
    n.breaks = 12,
    guide = guide_legend(reverse = TRUE)
  )+
  labs(
    fill = "Population\ncount",
    title = "Estimated population of Germany in 2020",
    subtitle = "Approx. 1km grid"
  )

Exercise 4_1: A simple map

Case study: Population dynamics

We will now dig deeper into population dynamics. We load a second layer which records the population size twenty years earlier - in 2000. We compare dimensions of both layers: they match.

ger_pop_2000 <- terra::rast("../../data/deu_ppp_2000_1km_Aggregated.tif") |> 
  terra::project("EPSG:25832")
  
ger_pop_2000
class       : SpatRaster 
dimensions  : 1183, 930, 1  (nrow, ncol, nlyr)
resolution  : 745.5165, 745.5165  (x, y)
extent      : 263554.6, 956885, 5235978, 6117924  (xmin, xmax, ymin, ymax)
coord. ref. : ETRS89 / UTM zone 32N (EPSG:25832) 
source(s)   : memory
name        : deu_ppp_2000_1km_Aggregated 
min value   :                       0.000 
max value   :                    3343.013 
ger_pop_2020
class       : SpatRaster 
dimensions  : 1183, 930, 1  (nrow, ncol, nlyr)
resolution  : 745.5165, 745.5165  (x, y)
extent      : 263554.6, 956885, 5235978, 6117924  (xmin, xmax, ymin, ymax)
coord. ref. : ETRS89 / UTM zone 32N (EPSG:25832) 
source(s)   : memory
name        : deu_ppp_2020_1km_Aggregated 
min value   :                       0.000 
max value   :                    4395.625

Case study: Population dynamics

We calculate simple differences for each cell between 2020 and 2000 to detect increases and decreases in estimated population per cell.

ger_pop_diff <- ger_pop_2020 - ger_pop_2000

terra::plot(ger_pop_diff)

Case study: Population dynamics

The state (Bundesland) of North Rhine-Westphalia looks quite interesting. The state has both growth and decline regions. We zoom into that region and explore dynamics within and between districts (Kreise).

NRW_districts <- sf::read_sf("../../data/VG250_KRS.shp") |> 
  filter(SN_L == "05")
  
ggplot(NRW_districts)+
  geom_sf(aes(fill = BEZ))

Case study: Population dynamics

The state (Bundesland) of North Rhine-Westphalia looks quite interesting. The state has both growth and decline regions. We zoom into that region and explore dynamics within and between districts (Kreise).

# Subset pop data to spatial extent of NRW
NRW_pop_diff <- terra::crop(
  ger_pop_diff,
  NRW_districts
) |> 
  terra::mask(NRW_districts)

# Let's visualize
ggplot()+
  geom_spatraster(data = NRW_pop_diff)+
  geom_sf(data = NRW_districts,
          fill = "transparent",
          color = "white",
          size = 2)

Case study: Population dynamics

We are now working with diverging values +/-0. The color palette should reflect that. We have several options to adjust it.

ggplot()+
  geom_spatraster(data = NRW_pop_diff)+
  geom_sf(data = NRW_districts,
          fill = "transparent",
          color = "white",
          size = 3)+
  scale_fill_viridis_c(option = "magma")

Case study: Population dynamics

We are now working with diverging values +/-0. The color palette should reflect that. We have several options to adjust it.

# Identify max and min to adjust for skewness of positive and negative values
max_val <- max(abs(minmax(NRW_pop_diff)))

ggplot()+
  geom_spatraster(data = NRW_pop_diff)+
  geom_sf(data = NRW_districts,
          fill = "transparent",
          color = "black",
          size = 3)+
  scale_fill_distiller(type = "div", 
                       palette = "PuOr", 
                       limits = c(-max_val, max_val)
                       )

Case study: Population dynamics

We are now working with diverging values +/-0. The color palette should reflect that. We have several options to adjust it.

# Or define manually for even greater control
ggplot()+
  geom_spatraster(data = NRW_pop_diff)+
  geom_sf(data = NRW_districts,
          fill = "transparent",
          color = "black",
          size = 3)+
  scale_fill_gradient2(
    low = "#2d004b",
    mid = "white",
    high = "#b35806",
    midpoint = 0
  )

Case study: Population dynamics

We also want to add text labels for the two districts with the highest population growth and decline.

# Identify growing and shrinking districts
NRW_districts <- NRW_pop_diff |> 
  extract(NRW_districts, fun = mean, na.rm =TRUE) |> 
  as_tibble() %>% 
  bind_cols(NRW_districts, .)

extremes <- NRW_districts |> 
  arrange(desc(deu_ppp_2020_1km_Aggregated)) |> 
  slice(c(1, n()))

Case study: Population dynamics

We also want to add text labels for the two districts with the highest population growth and decline.

ggplot()+
  geom_spatraster(data = NRW_pop_diff)+
  geom_sf(data = NRW_districts,
          fill = "transparent",
          color = "black",
          size = 3)+
  geom_sf_text(data = extremes, 
               aes(label = GEN), 
               size = 4, 
               color = "black"
               )+ 
  scale_fill_gradient2(
    low = "#2d004b",
    mid = "white",
    high = "#b35806",
    midpoint = 0
  )

Case study: Population dynamics

I prefer labels instead of text

library(ggrepel)

ggplot()+
  geom_spatraster(data = NRW_pop_diff)+
  geom_sf(data = NRW_districts,
          fill = "transparent",
          color = "black",
          size = 3)+
  geom_label_repel(
    data = extremes, 
    aes(geometry = geometry, label = GEN), 
    stat = "sf_coordinates",
    min.segment.length = 0, 
    segment.color = "black",   
    segment.size = 0.5,        
    size = 4,                  
    fill = "white",            
    color = "black",           
    label.size = 0.2           
    )+ 
  scale_fill_gradient2(
    low = "#2d004b",
    mid = "white",
    high = "#b35806",
    midpoint = 0
  )

Case study: Population dynamics

Let’s do some final polishing: 1. Turn NAs transparent, 2. add arrow and scale bar, 3. make final adjustments to the theme, and 4. adjust title, subtitle and caption.

library(ggspatial)

ggplot()+
  geom_spatraster(data = NRW_pop_diff)+
  geom_sf(data = NRW_districts,
          fill = "transparent",
          color = "black",
          size = 3)+
  geom_label_repel(
    data = extremes, 
    aes(geometry = geometry, label = GEN), 
    stat = "sf_coordinates",
    min.segment.length = 0, 
    segment.color = "black",   
    segment.size = 0.5,        
    size = 4,                  
    fill = "white",            
    color = "black",           
    label.size = 0.2           
  )+ 
  scale_fill_gradient2(
    low = "#2d004b",
    mid = "white",
    high = "#b35806",
    midpoint = 0,
    na.value = "transparent",
    name = "Population\nChange"
  )+
  theme_minimal()+
  theme(
    axis.title = element_blank()
  ) +
  labs(
    title = "Population change in NRW Districts (2000–2020)",
    subtitle = "In absolute numbers on 100x100m grid\nText labels identify highest growth and decline districts",
    caption = "Source: WorldPop (2018)"
  ) +
  annotation_scale(
    location = "br",
    width_hint = 0.3
  ) +
  annotation_north_arrow(
    location = "tl",
    which_north = "true",
    style = north_arrow_fancy_orienteering
  )

Beyond raw cell values: Contours

Let’s zoom into Cologne to explore the geom_spatraster() options a bit further.

cologne_bbox <- st_bbox(NRW_districts |> 
                          filter(GEN == "Köln")
                        )

# Subset pop data to spatial extent of Cologne
koeln_pop_2020 <- terra::crop(
  ger_pop_2020,
  cologne_bbox
)

ggplot()+
  geom_spatraster(data = koeln_pop_2020)+
  scale_fill_viridis_c(
    guide = guide_legend(reverse = TRUE)
  )+
  geom_sf(data = NRW_districts,
          fill = "transparent",
          color = "white",
          size = 8)+
  coord_sf(
    xlim = c(cologne_bbox["xmin"], 
             cologne_bbox["xmax"]),
    ylim = c(cologne_bbox["ymin"], 
             cologne_bbox["ymax"]),
    expand = FALSE
  )+
  theme_minimal()+
  theme(
    axis.title = element_blank()
  )

Beyond raw cell values: Contours

Let’s zoom into Cologne to explore the geom_spatraster options a bit further

# Plot contour
ggplot()+
  geom_spatraster_contour_filled(
    data = koeln_pop_2020
    )+
  scale_fill_viridis_d(
    guide = guide_legend(reverse = TRUE)
  )+
  geom_sf(data = NRW_districts,
          fill = "transparent",
          color = "white",
          size = 8)+
  coord_sf(
    xlim = c(cologne_bbox["xmin"], 
             cologne_bbox["xmax"]),
    ylim = c(cologne_bbox["ymin"], 
             cologne_bbox["ymax"]),
    expand = FALSE
  )+
  theme_minimal()+
  theme(
    axis.title = element_blank()
  )

Exercise 4_2: A fancy map